PACS: 71.15.-m, 71.27.+a, 75.10.Jm A Convergent Method for Calculating the Properties of Many Interacting Electrons

A method is presented for calculating binding energies and other properties of extended interacting systems using the projected density of transitions (PDoT) which is the probability distribution for transitions of different energies induced by a given localized operator, the operator on which the transitions are projected. It is shown that the transition contributing to the PDoT at each energy is the one which disturbs the system least, and so, by projecting on appropriate operators, the binding energies of equilibrium electronic states and the energies of their elementary excitations can be calculated. The PDoT may be expanded as a continued fraction by the recursion method, and as in other cases the continued fraction converges exponentially with the number of arithmetic operations, independent of the size of the system, in contrast to other numerical methods for which the number of operations increases with system size to maintain a given accuracy. These properties are illustrated with a calculation of the binding energies and zone-boundary spin-wave energies for an infinite spin-1/2 Heisenberg chain, which is compared with analytic results for this system and extrapolations from finite rings of spins. 1. Interacting Electrons Calculating binding energies for even a few interacting electrons is a formidable problem. Because of the interactions, the Schroedinger equation cannot be solved analytically, nor does it separate, and the number of parameters needed for a variational solution grows exponentially with the number of electrons. The ground state energy for a pair of interacting electrons which hop between just two spatial orbitals is already a substantial calculation. More electrons or orbitals requires a computer, and fewer than 100 interacting electrons defeat any numerical attempt to solve the Schroedinger equation variationally. The purpose of this paper is to show how properties of interacting electrons can be calculated so that the analytic or numerical effort required to obtain a given accuracy is independent of the number of electrons, and only grows logarithmically with the accuracy. This method arises from previous work on the states of non-interacting electrons in non-crystalline solids [1, 2, 3], and recent work [4] on the calculation of the projected density of transitions (PDoT) for interacting electrons. The main new result of this paper is that the binding energies and excitation energies of interacting electrons can be extracted from the PDoT, and this is illustrated with calculations for an infinite, spin-1/2 Heisenberg chain. While the Hamiltonian for a Heisenberg chain is simple, it still possesses all the complexity of interacting systems. For example the Hamiltonian for a chain of N+1 electrons has dimension 2N+1, one for each configuration of the spins, and so has 2N+1 stationary states whose energies per bond must lie between the ferromagnetic and anti-ferromagnetic limits. The number of states increases exponentially with the number of electrons while the energies of the states only increase linearly, so the density of states must also increase exponentially with the number of electrons. Even taking into account that the interaction conserves total spin, which can at most increase linearly with the number of electron, the density of states for a given total spin still increases exponentially with the number of electrons. In applying linear or non-linear variational methods to interacting systems, there must be a parameter for each degree of freedom which is to be varied independently, and so the number of variational parameters also increases exponentially with system size. Perturbative approaches diverge because the energy denominators decrease exponentially with the exponentially increasing number of coupled states in any energy interval. Other methods sample the states either statistically as in Monte Carlo, or by other criteria as in renormalization methods. The accuracy of Monte Carlo calculations only increases as the reciprocal of the root of the computational effort, and the errors in renormalization calculations depend on the criteria for neglecting states in a way which cannot be predicted in advance. The method presented here has its physics roots in the black body theorem of von Laue [5], which explains why the local (at a point) electromagnetic power spectrum is exponentially insensitive to its surroundings. Friedel [6] applied this idea to non-interacting electrons where the projected density of states (PDoS) the total density of state weighted by the probability of finding an electron in a particular orbital is exponentially insensitive to parts of the system distant from the orbital on which the states are projected. Friedel observed that moments of the PDoS could be calculated easily from powers of the electronic Hamiltonian, leading to local approaches to the electronic structure of solids [1] and to the recursion method for calculating projected densities of states from a continued fraction obtained by tridiagonalizing the electronic Hamiltonian. The focus on the PDoT and the use of energy-independent inner products, both of which are crucial for this approach, distinguish it from the projection methods developed by Zwanzig [7] and Mori [8] for correlation functions in classical and quantum statistical mechanics. While the stationary states of interacting systems have energies which grow with the size of the system, the stationary transitions, operators transforming one stationary state into another, need not have energies which depend on the system's size because they are the differences between energies of stationary states. The transitions satisfy Heisenberg's equation [4] which is equivalent to Schroedinger's equation, but does not become singular as the system gets larger, again because the transition energies do not increase. Just as non-interacting states are projected on a localized orbital for the PDoS, stationary transitions are projected on a localized operator for the PDoT which is then the total density of transitions weighted by the probability that each transition is induced by the localized operator on which it is being projected. It is crucial to the work presented here to avoid the singular behavior of stationary states with increasing system size. The mathematical roots of this work lie in the classical moment problem [9] and the properties of orthogonal polynomials [10]. The moment problem is that of reconstructing a distribution from its moments, integrals of the distribution over powers of the variable; in this work the distribution is the probability density for transitions and the variable is the energy of the transition. The surprising solution to the moment problem is that a continued fraction expansion of the distribution converges exponentially with the number of moments. Orthogonal polynomials in energy arise in this work as coefficients in expansions of the stationary transitions, and they are orthogonal with respect to integration over the PDoT. The numerical applications of these expansions are infinite dimensional analogs of the Lanczos method for finite matrices, in particular, convergence of the continued fraction expansion for the PDoS [11] corresponds to Paige's theorem [12] for the convergence of the Lanczos method. The rest of this paper is organized into 5 further Secs. In the next Sec. the local properties of the PDoT and its calculation by the recursion method are reviewed briefly. In addition to the review, there is a construction of density matrices for stationary transitions from which expectation values can be calculated. In Sec. 3, the property of the PDoT that the transition which contributes at each energy is the one which is qualitatively most localized is used to show how binding energies, excitation energies, and other quantities can be calculated for individual equilibrium states. The way transitions with different localization properties contribute to the PDoT is illustrated in Sec. 4 which contains analytic examples of a localized transition which is degenerate with a band of delocalized transitions, and an example of two degenerate bands of transitions with different localization properties. In Sec. 5 is the description of a calculation of the equilibrium binding energies and zone boundary spin-wave energies for electrons in an infinite Heisenberg chain, which is compared to analytic results for this system and some numerical results for finite rings. The last Sec. contains a discussion of how this method is related to thermodynamical and statistical mechanical approaches. While the Heisenberg chain is used throughout this paper to illustrate the method presented, it is intended to be clear that this method applies to a wide range of systems. 2. Projected Density of Transitions, the Recursion Method, and Expectation Values This work is based on the two local properties of the PDoT analogous to properties of the PDoS: that it is exponentially insensitive to distant parts of the interacting system, and that the transition which contributes to the PDoT at each energy is the one which is most localized, and so changes the system least. Because of this, there follows a brief review of the relation of the PDoT to the energy resolvent for the evolution of operators, its calculation by the recursion method, and an expansion of the density matrices for stationary transitions. For a system with N non-degenerate (to avoid complications) stationary energies {Eα} and corresponding states {ψα(r)}, expressed as functions of r which stands for the dynamical variables such as position, spin, and so on. If u is an operator on which the transitions are to be projected, for example the annihilation operator for an electron in a particular orbital, then the PDoT for this system is defined to be, ( ) ( ) α β β α α β + δ ∑ ∫ ψ ψ = E E E d u * ) ( N) (1/ (E) g 2